home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / debug_fft.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  9.5 KB  |  438 lines
  1. #define MAX_SIZE 1111
  2.  
  3. *****************************************************************************
  4. *                                        *
  5. *    1D (slow) Fourier Transform                        *
  6. *                                        *
  7. *****************************************************************************
  8.       subroutine dft1du( plusminus, n, a, inc)
  9.       double precision a(*)
  10.       double precision w(MAX_SIZE)
  11.       integer n, plusminus, inc
  12.       
  13.       double precision u, v, tpi
  14.       TPI = 2 * 3.14159265358979323846
  15.  
  16.       if( inc .ne. 1) goto 11
  17.       if ( plusminus .eq. -1) then
  18.     l = (n+2)/2
  19. c$doacross local(k,u,v)
  20.     do k = 1, l
  21.         u = 0.0
  22.         v = 0.0
  23.         do i = 1, n
  24.         u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
  25.         v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
  26.         end do
  27.         w((2*k-2)+1) = u
  28.         w((2*k-1)+1) = v
  29.     end do
  30. c$doacross
  31.     do i = 1, 2*l
  32.         a(i) = w(i)
  33.     end do
  34.       elseif (plusminus .eq. 1) then
  35.     Print *,'NOT Implemented'
  36.       else
  37.     print *, 'Error first argument of DFT1DU should -1 or +1'
  38.     stop
  39.  11    print *, 'This Simple version support only INC = 1'
  40.     stop
  41.       endif
  42.  
  43.       return
  44.       end
  45.  
  46.  
  47.       subroutine zft1d( plusminus, n, a, inc)
  48.       implicit double precision ( a-h, o-z)
  49.       double complex a(*)
  50.       double complex w(MAX_SIZE)
  51.       integer n, plusminus
  52.       
  53.       double complex cp
  54.       double precision re, im, tpi
  55.  
  56.       tpi = 2.*3.14159265358979323846
  57.  
  58.     do i = 1, n
  59.         w(i) = 0.0
  60.         do k = 1, n
  61.         im = mod(((i-1)*(k-1)), n)
  62.         im = dfloat( plusminus*im)*tpi / dfloat(n) 
  63.         re = cos( im)
  64.         im = sin( im)
  65.         cp = dcmplx( re, im)
  66.         w(i) = w(i) + cp*a(inc*(k-1)+1)
  67.         end do
  68.     end do
  69.     do i = 1, n
  70.         a(inc*(i-1)+1) = w(i)
  71.     end do
  72.  
  73.       return
  74.       end
  75.  
  76.  
  77.  
  78. *****************************************************************************
  79. *                                        *
  80. *    2D (slow) Fourier Transform                        *
  81. *                                        *
  82. *****************************************************************************
  83.       subroutine zft2d( plusminus, n1, n2, w, ldw)
  84.       double complex w(ldw,n2)
  85.       integer n1, n2, lda, ldw, plusminus
  86.  
  87. c$doacross local(j)
  88.       do j = 1, n2
  89.     call zft1d( plusminus, n1, w(1,j), 1)
  90.       end do
  91.  
  92. c$doacross local(i)
  93.       do i = 1, n1
  94.     call zft1d( plusminus, n2, w(i,1), ldw)
  95.       end do
  96.  
  97.       return
  98.       end
  99.  
  100.       subroutine dft2du( plusminus, n1, n2, w, ldw)
  101.       Double Precision w(ldw,n2)
  102.       integer n1, n2, lda, ldw, plusminus, l, lh
  103.  
  104.       l = ((n1+2)/2)*2
  105.       lh = ldw/2
  106.       if (ldw .lt. l) then
  107.     PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
  108.       endif
  109.       if( mod(ldw,2) .ne. 0) then
  110.     PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
  111.       endif 
  112.       if (plusminus .eq. -1) then
  113. c$doacross local(j)
  114.     do j = 1, n2
  115.         call dft1du( plusminus, n1, w(1,j), 1)
  116.     end do
  117. c$doacross local(i,j,my)
  118.     do i = 1, l, 2
  119.         call zft1d( plusminus, n2, w(i,1), lh)
  120.     end do
  121.       else
  122. c$doacross local(i,j,my)
  123.     do i = 1, l, 2
  124.         call zft1d( plusminus, n2, w(i,1), lh)
  125.     end do
  126. c$doacross local(j)
  127.     do j = 1, n2
  128.         call dft1du( plusminus, n1, w(1,j), 1)
  129.     end do
  130.       end if
  131.       return
  132.       end
  133.       
  134.  
  135. *****************************************************************************
  136. *
  137. *    3D (slow) Fourier Transform
  138. *
  139. *****************************************************************************
  140.  
  141.       subroutine zft3d(plusminus,n1,n2,n3,w,ld1,ld2)
  142.       double complex w(ld1,ld2,n3)
  143.       integer n1,n2,n3,ld1,ld2,plusminus
  144. c
  145. c   transform along X  first ...
  146. c
  147. c$doacross local(i,j,k)
  148.       do k = 1, n3
  149.        do j = 1, n2
  150.         call zft1d( plusminus, n1, w(1,j,k), 1)
  151.        end do
  152.       end do
  153. c
  154. c   transform along Y then ...
  155. c
  156. c$doacross local(i,j,k)
  157.       do k = 1,n3
  158.        do i = 1,n1
  159.         call zft1d( plusminus, n2, w(i,1,k), ld1)
  160.        end do
  161.       end do
  162. c
  163. c   transform along Z finally ...
  164. c
  165. c$doacross local(i,j,k)
  166.       do i = 1, n1
  167.        do j = 1, n2
  168.         call zft1d( plusminus, n3, w(i,j,1), ld1*ld2)
  169.        end do
  170.       end do
  171.  
  172.       return
  173.       end
  174.  
  175.       subroutine dft3du(plusminus,n1,n2,n3,w,ld1,ld2)
  176.       Double Precision w(ld1,ld2,n3)
  177.       integer n1,n2,n3,ld1,ld2,plusminus, l
  178.  
  179.       if( plusminus .ne. -1) then
  180.     print *, 'Option non implemented YET'
  181.       end if
  182.       if( mod(ld1,2) .ne. 0) then
  183.     PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
  184.       endif 
  185.       l = ((n1+2)/2)*2
  186.       if( ld1 .lt. l ) then
  187.     print *, 'Leading dimension not sufficient'
  188.       endif
  189. c
  190. c   trandform along X  first ...
  191. c
  192. c$doacross local(i,j,k)
  193.       do k = 1, n3
  194.        do j = 1, n2
  195.         call dft1du( plusminus, n1, w(1,j,k), 1)
  196.        end do
  197.       end do
  198. c
  199. c   trandform along Y then ...
  200. c
  201. c$doacross local(i,k)
  202.       do k = 1,n3
  203.     do i = 1, l, 2
  204.         call zft1d( plusminus, n2, w(i,1,k), ld1/2)
  205.     end do
  206.       end do
  207. c
  208. c   trandform along Z finally ...
  209. c
  210. c$doacross local(i,j)
  211.       do j = 1,n2
  212.     do i = 1, l, 2
  213.         call zft1d( plusminus, n3, w(i,j,1), (ld1*ld2)/2)
  214.     end do
  215.       end do
  216.  
  217.       return
  218.       end
  219.  
  220. *****************************************************************************
  221. *                                        *
  222. *    1D (slow) Fourier Transform                        *
  223. *                                        *
  224. *****************************************************************************
  225.       subroutine sft1du( plusminus, n, a, inc)
  226.       real a(*)
  227.       real w(MAX_SIZE)
  228.       integer n, plusminus, inc
  229.       
  230.       real u, v, tpi
  231.       TPI = 2 * 3.14159265358979323846
  232.  
  233.       if( inc .ne . 1) goto 11
  234.       if ( plusminus .eq. -1) then
  235.     l = (n+2)/2
  236. c$doacross local(k,u,v)
  237.     do k = 1, l
  238.         u = 0.0
  239.         v = 0.0
  240.         do i = 1, n
  241.         u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
  242.         v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
  243.         end do
  244.         w((2*k-2)+1) = u
  245.         w((2*k-1)+1) = v
  246.     end do
  247. c$doacross
  248.     do i = 1, 2*l
  249.         a(i) = w(i)
  250.     end do
  251.       elseif (plusminus .eq. 1) then
  252.     print *,'NOT Implemented'
  253.       else
  254.     print *, 'Error first argument of DFT1DU should -1 or +1'
  255.     stop
  256.  11    print *, 'This Simple version support only INC = 1'
  257.     stop
  258.       endif
  259.  
  260.       return
  261.       end
  262.  
  263.  
  264.       subroutine cft1d( plusminus, n, a, inc)
  265.       implicit real ( a-h, o-z)
  266.       complex a(*)
  267.       complex w(MAX_SIZE)
  268.       integer n, plusminus
  269.       
  270.       complex cp
  271.       real re, im, tpi
  272.  
  273.       tpi = 2.*3.14159265358979323846
  274.  
  275.     do i = 1, n
  276.         w(i) = 0.0
  277.         do k = 1, n
  278.         im = mod(((i-1)*(k-1)), n)
  279.         im = float( plusminus*im)*tpi / float(n) 
  280.         re = cos( im)
  281.         im = sin( im)
  282.         cp = cmplx( re, im)
  283.         w(i) = w(i) + cp*a(inc*(k-1)+1)
  284.         end do
  285.     end do
  286.     do i = 1, n
  287.         a(inc*(i-1)+1) = w(i)
  288.     end do
  289.  
  290.       return
  291.       end
  292.  
  293.  
  294.  
  295. *****************************************************************************
  296. *                                        *
  297. *    2D (slow) Fourier Transform                        *
  298. *                                        *
  299. *****************************************************************************
  300.       subroutine cft2d( plusminus, n1, n2, w, ldw)
  301.       complex w(ldw,n2)
  302.       integer n1, n2, lda, ldw, plusminus
  303.  
  304. c$doacross local(j)
  305.       do j = 1, n2
  306.     call cft1d( plusminus, n1, w(1,j), 1)
  307.       end do
  308.  
  309. c$doacross 
  310.       do i = 1, n1
  311.     call cft1d( plusminus, n2, w(i,1), ldw)
  312.       end do
  313.  
  314.       return
  315.       end
  316.  
  317.       subroutine sft2du( plusminus, n1, n2, w, ldw)
  318.       real w(ldw,n2)
  319.       integer n1, n2, lda, ldw, plusminus, l, lh
  320.  
  321.       l = ((n1+2)/2)*2
  322.       lh = ldw/2
  323.       if (ldw .lt. l) then
  324.     PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
  325.       endif
  326.       if( mod(ldw,2) .ne. 0) then
  327.     PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
  328.       endif 
  329.       if (plusminus .eq. -1) then
  330. c$doacross local(j)
  331.     do j = 1, n2
  332.         call sft1du( plusminus, n1, w(1,j), 1)
  333.     end do
  334.  
  335. c$doacross local(i,j,my)
  336.     do i = 1, l, 2
  337.         call cft1d( plusminus, n2, w(i,1), lh)
  338.     end do
  339.       else
  340. c$doacross local(i,j,my)
  341.     do i = 1, l, 2
  342.         call cft1d( plusminus, n2, w(i,1), lh)
  343.     end do
  344. c$doacross local(j)
  345.     do j = 1, n2
  346.         call sft1du( plusminus, n1, w(1,j), 1)
  347.     end do
  348.       end if
  349.       return
  350.       end
  351.       
  352.  
  353. *****************************************************************************
  354. *
  355. *    3D (slow) Fourier Transform
  356. *
  357. *****************************************************************************
  358.  
  359.       subroutine cft3d(plusminus,n1,n2,n3,w,ld1,ld2)
  360.       complex w(ld1,ld2,n3)
  361.       integer n1,n2,n3,ld1,ld2,plusminus
  362. c
  363. c   transform along X  first ...
  364. c
  365. c$doacross local(i,j,k)
  366.       do k = 1, n3
  367.        do j = 1, n2
  368.         call cft1d( plusminus, n1, w(1,j,k), 1)
  369.        end do
  370.       end do
  371. c
  372. c   transform along Y then ...
  373. c
  374. c$doacross local(i,j,k)
  375.       do k = 1,n3
  376.        do i = 1,n1
  377.         call cft1d( plusminus, n2, w(i,1,k), ld1)
  378.        end do
  379.       end do
  380. c
  381. c   transform along Z finally ...
  382. c
  383. c$doacross local(i,j,k)
  384.       do i = 1, n1
  385.        do j = 1, n2
  386.         call cft1d( plusminus, n3, w(i,j,1), ld1*ld2)
  387.        end do
  388.       end do
  389.  
  390.       return
  391.       end
  392.  
  393.       subroutine sft3du(plusminus,n1,n2,n3,w,ld1,ld2)
  394.       real w(ld1,ld2,n3)
  395.       integer n1,n2,n3,ld1,ld2,plusminus, l
  396.  
  397.       if( plusminus .ne. -1) then
  398.     print *, 'Option non implemented YET'
  399.       end if
  400.       if( mod(ld1,2) .ne. 0) then
  401.     PRINT *, 'Warning ... this Checking routine assumes an even leading Dimension '
  402.       endif 
  403.       l = ((n1+2)/2)*2
  404.       if( ld1 .lt. l ) then
  405.     print *, 'Leading dimension not sufficient'
  406.       endif
  407. c
  408. c   trandform along X  first ...
  409. c
  410. c$doacross local(i,j,k)
  411.       do k = 1, n3
  412.        do j = 1, n2
  413.         call sft1du( plusminus, n1, w(1,j,k), 1)
  414.        end do
  415.       end do
  416. c
  417. c   transform along Y then ...
  418. c
  419. c$doacross local(i,k)
  420.       do k = 1,n3
  421.     do i = 1, l, 2
  422.         call cft1d( plusminus, n2, w(i,1,k), ld1/2)
  423.     end do
  424.       end do
  425. c
  426. c   transform along Z finally ...
  427. c
  428. c$doacross local(i,j)
  429.       do j = 1,n2
  430.     do i = 1, l, 2
  431.         call cft1d( plusminus, n3, w(i,j,1), (ld1*ld2)/2)
  432.     end do
  433.       end do
  434.  
  435.       return
  436.       end
  437.  
  438.